import struct, os
import dbf

# ESRI shapefile stores integer and double-precision numbers
# Integer = signed 32-bit integer (4 bytes)
# Double = signed 64-bit IEEE double-precision floating point number (8 bytes)
# NoData = floating point smaller than -10e38

# the bytes in the file header and record contents in main file
# are in little-endian byte order. everything else is big-endian.

# main file header is 100 bytes long.

# Position  Field          Value        Type      Byte Order
#-------------------------------------------------------------
# Byte 0    File Code      9994         Integer   Big   
# Byte 4    Unused         0            Integer   Big   
# Byte 8    Unused         0            Integer   Big   
# Byte 12   Unused         0            Integer   Big   
# Byte 16   Unused         0            Integer   Big   
# Byte 20   Unused         0            Integer   Big   
# Byte 24   File Length    File Length  Integer   Big   
# Byte 28   Version        1000         Integer   Little
# Byte 32   Shape Type     Shape Type   Integer   Little
# Byte 36   Bounding Box   Xmin         Double    Little
# Byte 44   Bounding Box   Ymin         Double    Little
# Byte 52   Bounding Box   Xmax         Double    Little
# Byte 60   Bounding Box   Ymax         Double    Little
# Byte 68*  Bounding Box   Zmin         Double    Little
# Byte 76*  Bounding Box   Zmax         Double    Little
# Byte 84*  Bounding Box   Mmin         Double    Little
# Byte 92*  Bounding Box   Mmax         Double    Little
            
# * Unused, with value 0.0, if not Measured or Z type

def shapetype_to_string(shapetype):
    names = {0:'Null Shape', 1:'Point', 3:'PolyLine', 5:'Polygon',
             8:'MultiPoint', 11:'PointZ', 13:'PolyLineZ', 15:'PolygonZ',
             18:'MultiPointZ', 21:'PointM', 23:'PolyLineM', 25:'PolygonM',
             28:'MultiPointM', 31:'MultiPatch'}
    if shapetype in names:
        return names[shapetype]
    return None
    

def unpack_bounding_box(s):
    return struct.unpack('<4d', s)

def unpack_point(s):
    return struct.unpack('<2d', s)


class BoundingBox:
    __slots__ = ['xmin','ymin','xmax','ymax']

class Point:
    __slots__ = ['id', 'x','y']
    def __init__(self, id, x, y):
        self.id = id
        self.x = x
        self.y = y

class Polyline:
    __slots__ = ['id', 'vertices']
    def __init__(self, id, vertices):
        self.id = id
        self.vertices = vertices

class Polygon:
    __slots__ = ['id', 'vertices']
    def __init__(self, id, vertices):
        self.id = id
        self.vertices = vertices


class Shapefile(object):
    def __init__(self, filename):
        self.filename = filename

        # stored geometry
        self.points = []
        self.polylines = []
        self.polygons = []

        with open(filename, 'rb') as f:
            self.read_file_header(f)
            print "filename:", filename
            print "shape type:", shapetype_to_string(self.shape_type)
            print "bounds:", self.xmin, self.ymin, self.xmax, self.ymax

            index = 0
            while self.read_record(index, f):
                index += 1

        # load the attributes from the .dbf file
        #root, ext = os.path.splitext(filename)
        #self.attribute_table = dbf.dBaseTable(root + '.dbf')
        

    def read_file_header(self, f):
        headerbytes = f.read(100)
        if headerbytes:
            self.file_code, = struct.unpack('>i', headerbytes[0:4])
            assert self.file_code==9994
            self.file_length, = struct.unpack('>i', headerbytes[24:28])
            version, shape_type, xmin, ymin, xmax, ymax, zmin, zmax, mmin, mmax = struct.unpack('<2i8d', headerbytes[28:])
            self.version = version
            assert self.version==1000
            self.shape_type = shape_type
            self.xmin = xmin
            self.xmax = xmax
            self.ymin = ymin
            self.ymax = ymax
            self.zmin = zmin
            self.zmax = zmax
            self.mmin = mmin
            self.mmax = mmax

    def read_record(self, objectid, f):
        record_header = f.read(8)
        if record_header:
            record_number, length = struct.unpack('>2i', record_header)
        else:
            return None
        record = f.read(length*2)
        shape_type, = struct.unpack('<i', record[0:4])
        #print "Record:", record_number, length, shapetype_to_string(shape_type)
        if shape_type == 0: # null shape
            return None

        elif shape_type == 1: # point
            x, y = struct.unpack('<2d', record[4:20])
            point = Point(objectid, x, y)
            self.points.append(point)
            return point
        
        elif shape_type == 8: # multipoint
            xmin, ymin, xmax, ymax = struct.unpack('<4d', record[4:36])
            num_points, = struct.unpack('<i', record[36:40])
            point_list = []
            pos = 40
            for i in range(num_points):
                x, y = struct.unpack('<2d', record[pos:pos+16])
                point_list.append(Point(objectid, x, y))
                pos += 16
            self.points.extend(point_list)
            return point_list

        elif shape_type == 3: # polyline
            xmin, ymin, xmax, ymax = struct.unpack('<4d', record[4:36])
            num_parts, = struct.unpack('<i', record[36:40])
            num_points, = struct.unpack('<i', record[40:44])
            parts = struct.unpack('<'+str(num_parts)+'i', record[44:44+num_parts*4])
            point_list = []
            pos = 44 + num_parts*4
            for i in range(num_points):
                point = struct.unpack('<2d', record[pos:pos+16])
                point_list.append(point)
                pos += 16
            # separate point_list into polyline parts
            polylines = []
            parts = parts + (num_points,)
            for i in range(num_parts):
                start = parts[i]
                stop = parts[i+1]
                polyline = Polyline(objectid, point_list[start:stop])
                polylines.append( polyline )
                #print 'PolyLine %d, %d:' % (num_parts, num_points), polyline
            self.polylines.extend(polylines)
            return polylines

        elif shape_type == 5: # polygon
            # Polygons are made up of parts called "rings".
            # Rings are sequences of vertices.  The first and last
            # vertex must be the same, even though this is redundant.
            # The interior of the polygon is on the right edge.
            # So clockwise ordering defines interior and
            # counterclockwise is used for defining holes.
            # Want to use a constrained delauney triangulation to break
            # the concave polygons into triangles.
            xmin, ymin, xmax, ymax = struct.unpack('<4d', record[4:36])
            num_parts, = struct.unpack('<i', record[36:40])
            num_points, = struct.unpack('<i', record[40:44])
            parts = struct.unpack('<'+str(num_parts)+'i', record[44:44+num_parts*4])
            point_list = []
            pos = 44 + num_parts*4
            for i in range(num_points):
                point = struct.unpack('<2d', record[pos:pos+16])
                point_list.append(point)
                pos += 16
            # separate point_list into polygon parts
            polygons = []
            parts = parts + (num_points,)
            for i in range(num_parts):
                start = parts[i]
                stop = parts[i+1]
                # Don't include the last vertex, since same as the first.
                polygon = Polygon(objectid, point_list[start:stop-1])
                polygons.append( polygon )
                #print 'Polygon %d, %d:' % (num_parts, num_points), polyline
            self.polygons.extend(polygons)
            return polygons
            
